// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#include "transportio.h"

using namespace std;

// 
// lead struct contains everything that is the property of the lead. 
//
// The input parameters that define the lead are : 'tol' 'h'  's'  'v'  't'
// 'tol'is the channel direction and open-close separation tolerance (~1e-8)
// 'h' is the hamiltonian of one unit cell and 's' is the corresponding overlap matrix.
// 'v' is the coupling hamiltonian between two neighbouring unit cell and 't' is the corrsponding overlap matrix.
//
// The output parameters that are defined by the lead are : 'sl'  'sr'  'wlo'  'iwlo'  'wro'  'iwro'
// where l (r) refers to the left (right) 
// 'sl' 'sr' the left and right self energy [square matrix]
// 'wlo' 'wro' the left and right moving eigen solutions  [not square matrix in general]
// 'iwlo' 'iwro' the corresponding left duals i.e. iwlo*wlo and iwro*wro gives unit matrix 
// 
struct lead
{
    //input
    double tol;
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> h,v,s,t;    
    //computational
    int n,noocl,noocr,nocl,nocr;
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> id,zero,he,ve,ive,ave,w,z,wl ,wr,iwl,iwr,wvl,wvr,iwvl,iwvr,dl,idl,dr,idr,d;
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,1> gv,gvl,gvr,zl,zr;  
    Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> gvml,gvmr,igvml,igvmr;
    Eigen::Matrix<int,Eigen::Dynamic,1> index;
    //output
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> sl,sr,wlo,iwlo,wro,iwro;
    //methods
    bool self_energy(double e);
    //label
    std::string label;
    //option
    bool advanced;
    
};

bool lead::self_energy(double e)
{
	std::complex<double> ii(0,1.0);           // imaginary constant
	
	n = h.rows();                             // n is the dimension of the unit cell hamiltonian
		 
	id   = Eigen::MatrixXcd::Identity(n,n);    // indentity matrix with dimension n
	zero = Eigen::MatrixXcd::Zero(n,n);        // zero matrix with dimension n
	
	he = s*e-h;                               // taking into account the overlaps
	ve = t*e-v;                               // ... for the coupling hamiltonian

#ifdef DEBUG
        std::cout << "energy: " << e << std::endl;


        std::cout << "h:" << std::endl;
	out_complex_matrix(h,std::cout);
        std::cout << "s:" << std::endl;
	out_complex_matrix(s,std::cout);

        std::cout << "v:" << std::endl;
	out_complex_matrix(v,std::cout);
        std::cout << "t:" << std::endl;
	out_complex_matrix(t,std::cout);



        std::cout << "he:" << std::endl;
	out_complex_matrix(he,std::cout);
        std::cout << "ve:" << std::endl;
	out_complex_matrix(ve,std::cout);
#endif


	//ive = ve.inverse();                     // inverse of the e-dependent coupling [assumed is possible, if not then forget this code! or use randomize]
        Eigen::FullPivLU<Eigen::MatrixXcd> lu(ve);
        if( lu.isInvertible() ) ive = lu.inverse();
        else throw std::runtime_error("Coupling matrix is singular. Use randomize options!");
           
           
                   
	ave = ve.adjoint();                       // adjoint ov the e-dependent coupling

#ifdef DEBUG
        std::cout << "ave:" << std::endl;
	out_complex_matrix(ave,std::cout);
        std::cout << "ive:" << std::endl;
	out_complex_matrix(ive,std::cout);
#endif


	d.resize(2*n,2*n);                        // constructing the 2n x 2n matrix for second order eigenvalue problem
	d.block(0,0,n,n) = -(ive*he);      // ...
	d.block(0,n,n,n) = -(ive*ave);     // ... 
	d.block(n,0,n,n) = id;             // ...
	d.block(n,n,n,n) = zero;           // ...

#ifdef DEBUG
        std::cout << "d:" << std::endl;
	out_complex_matrix(d,std::cout);
#endif

        Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(d);  // solving eigenvalue problem
	w = eigen_solver.eigenvectors().block(0,0,n,2*n);                                              // 2n column eigenvectors with n rows; 
	z = eigen_solver.eigenvalues();                                                                // 2n eigenvalues
	
        for(int i = 0; i < 2*n; i++)              // normalizing the 2n eigenvectors
		w.col(i).normalize();             // ...
        for(int i = 0; i < 2*n; i++)              // normalizing the 2n eigenvectors
	{
                std::complex<double> overlap=w.col(i).dot((t.adjoint()/z(i)+s+t*z(i))*w.col(i));
        	w.col(i) /= sqrt(overlap);        // ...
        }

#ifdef DEBUG
        std::cout << "w:" << std::endl;
	out_complex_matrix(w,std::cout);
        std::cout << "z:" << std::endl;
	out_complex_matrix(z,std::cout);
#endif
	gv.resize(2*n);                           // calculating group velocities
	for(int i = 0; i<2*n; i++)                // ...
	gv(i) = -ii*w.col(i).dot((ave/z(i)-ve*z(i))*w.col(i));   // groupvelocity formula

/*
	for(int i = 0; i<2*n; i++)                // ...
	{
		std::complex<double> overlap=w.col(i).dot((t.adjoint()/z(i)+s+t*z(i))*w.col(i));       // overlap contributions
		//gv(i) = -ii*w.col(i).dot((ave/z(i)-ve*z(i))*w.col(i))/overlap;                         // groupvelocity formula
		gv(i) = -ii*w.col(i).dot((ave/z(i)-ve*z(i))*w.col(i));                         // groupvelocity formula
	}
*/

#ifdef DEBUG
	std::cout << "group velocity:" << std::endl;
	out_complex_vector(gv,std::cout);
#endif

	if (advanced) out_complex_vector_as_line(label+"-gv",e,gv,std::cout);


	wl.resize(n,n);                           // only left going eigen solutions;  each column a vector  
	wr.resize(n,n);                           // only right going solutions
	gvl.resize(n);                            // group velcoties to left
	gvr.resize(n);                            // group velocities to right
	zl.resize(n);                             // eigenvalues for left going solutions
	zr.resize(n);                             // eigenvalues for right going solutions

	//double tol = 0.00001;                     // channel direction and open-close separation tolerance
	int kl = 0; int kr = 0;                   // selecting apart the left and right going quantities
	for(int i = 0; i < 2*n; i++)                                                                 if( ( abs(abs(z(i))-1.0) <= tol) && (real(gv(i)) < 0.0) )     
        { for(int j = 0; j < n; j++) wl(j,kl) = w(j,i); zl(kl) = z(i); gvl(kl) = gv(i); kl++; } else if( ( abs(abs(z(i))-1.0) <= tol) && (real(gv(i)) > 0.0) )
        { for(int j = 0; j < n; j++) wr(j,kr) = w(j,i);	zr(kr) = z(i); gvr(kr) = gv(i); kr++; }
	noocl = kl; noocr = kr;                   // number of left/right open channels; they should be the same
	for(int i = 0; i <2*n; i++)                                                                  if( ( abs(abs(z(i))-1.0) >  tol) && (abs(z(i)) > 1.0) )   
        { for(int j = 0; j < n; j++) wl(j,kl) = w(j,i); zl(kl) = z(i); gvl(kl) = gv(i); kl++; } else if( ( abs(abs(z(i))-1.0) >  tol) && (abs(z(i)) < 1.0) ) 
        { for(int j = 0; j < n; j++) wr(j,kr) = w(j,i);	zr(kr) = z(i); gvr(kr) = gv(i); kr++; }
	nocl = kl; nocr = kr;                     // number of left/right channels (all);  they should be the same
        if( (noocl != noocr) || (nocl != nocr) ) throw "Channel partitioning error! (Check transport input, adjust tolerance or try randomize.)";
        if(noocl*noocr == 0) return false;

#ifdef DEBUG
	std::cout << "wl:" << std::endl;
        out_complex_matrix(wl,std::cout);
	std::cout << "wr:" << std::endl;
        out_complex_matrix(wr,std::cout);

	std::cout << "back subtitution to the original bloch equation:" << std::endl;
	out_complex_matrix(ave*wl*zl.asDiagonal().inverse()*wl.inverse() + he + ve*wl*zl.asDiagonal()*wl.inverse(),std::cout);
	out_complex_matrix(ave*wr*zr.asDiagonal().inverse()*wr.inverse() + he + ve*wr*zr.asDiagonal()*wr.inverse(),std::cout);
#endif
	iwl = wl.inverse();                       // calculating left duals; n row vectors
	iwr = wr.inverse();                       // calculating right duals; n row vectors

	//dl  = wl*zl.asDiagonal()*iwl;           // calculating the two solutions of the second order matrix equations
	idl = wl*zl.asDiagonal().inverse()*iwl;   // ...
	dr  = wr*zr.asDiagonal()*iwr;             // ...
	//idr = wr*zr.asDiagonal().inverse()*iwr; // ...

	sl=-ave*idl;                              // calculating left self energy
	sr=-ve*dr;                                // ... and right self energy

#ifdef DEBUG
        std::cout << "back substituting into the recursion equations:" << std::endl;
	out_complex_matrix(ave*(he-sl).inverse()*ve-sl,std::cout);
	out_complex_matrix(ve*(he-sr).inverse()*ave-sr,std::cout);
        std::cout << "reduced hamiltonian (with self-energies)" << std::endl;
	out_complex_matrix(he -sl - sr,std::cout);
#endif
	igvml=gvl.array().abs().sqrt().inverse(); // calculating (inverse) groupvelocity; for unit-flux normalization
	igvmr=gvr.array().abs().sqrt().inverse(); // ...
	gvml=gvl.array().abs().sqrt();            // ...
	gvmr=gvr.array().abs().sqrt();            // ...
	
	wvl = wl*igvml.asDiagonal();              // unit-flux normalization
	wvr = wr*igvmr.asDiagonal();              // ...
	
        iwvl = gvml.asDiagonal()*iwl;             // unit-flux normalization for the duals
	iwvr = gvmr.asDiagonal()*iwr;             // ...
	
	wlo = wvl.block(0,0,n,noocl);             // keeping only the left open channels
	iwlo = iwvl.block(0,0,noocl,n);           // ... for the duals as well
	
	wro = wvr.block(0,0,n,noocr);             // keeping only the right open channels 
	iwro = iwvr.block(0,0,noocr,n);           // ... for the duals as well

#ifdef DEBUG
	std::cout << "checking duals for open channels:" << std::endl;
	out_complex_matrix(iwlo*wlo,std::cout);
        out_complex_matrix(iwro*wro,std::cout);
#endif
        return true;
}

//
// scatterer struct contains everything that is the property of the scatterer and its interaction with leads
//
// The input parameters that define the scatterer are : 'h' 's' 'va' 'ta' 'vb' 'tb'
// 'h' is the hamiltonian of the scatterer and 's' is the corresponding overlap matrix.
// 'va' is the coupling hamiltonian between the a-lead and the scatterer and 'ta' is the corresponding overlap matrix.
// 'vb' is the coupling hamiltonian between the b-lead and the scatterer and 'tb' is the corresponding overlap matrix.
//
// The output parameters that are defined by the scatterer are : 'sa' 'sb' 'sab'
// 'sa' and 'sb' the self energy towards a-lead and b-lead.
// 'sab' "self" energy like term between the two leads
//
struct scatterer
{
    //input
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> h,s,va,ta,vb,tb; 
    //computational
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> he,vae,vbe,ihe; //,ihvae,ihvbe; 
    //output
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> sa,sb,sab; 
    //methods
    void self_energy(double e);
};

void scatterer::self_energy(double e)
{
	std::complex<double> ii(0,1.0);           // imaginary constant

	he = s*e-h;                               // taking into account the overlaps
	vae = ta*e-va;                            // ... for the coupling hamiltonians
        vbe = tb*e-vb;                            // ...
        
        //he += s*eps*ii;                           // imaginary 
        ihe = he.inverse();                       // inverse of e-dependent hamiltonian
/*
        ihvae = ihe*vae.adjoint();                // matrices to project lead wave funcs
        ihvbe = ihe*vbe;                          // ...  to scatterer wave funcs.
*/               
        sa = vae*ihe*vae.adjoint();               // self energy towards a-lead
        sb = vbe.adjoint()*ihe*vbe;               // ... b-lead
        sab = vae*ihe*vbe;                        // ... between a-lead and b-lead
}

struct scattering
{
    //computational
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> sm,ig,dsm,g,kaa,kbb,kab,kba,k,id,ttab;
    //output
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> raa,rbb,tab,tba,s,xaar,xabl,xbar,xbbl,xmr,xml,yml,ymr,ttabev;
    Eigen::VectorXcd ttabe;
    double coefaa,coefbb,coefab,coefba,checknorm,maxcoefab,maxcoefba;
    int ano,bno;
    Eigen::Matrix<std::complex<double>,Eigen::Dynamic,1> agvl,agvr,bgvl,bgvr;
    //methods
    bool transport(double e,lead & a,scatterer & m,lead & b);
};

bool scattering::transport(double e,lead & a,scatterer & m,lead & b)
{
        if(a.self_energy(e)) ano=a.noocl; else ano=0;              // no transport if no open channel
        if(b.self_energy(e)) bno=b.noocl; else bno=0;              // ...
        if(ano*bno == 0) return false;                             // ...

        m.self_energy(e);                         // calculating scatterer properties
        
        sm.resize(a.n+b.n,a.n+b.n);               // scatterer selfenergy
        sm.block(0,0,a.n,a.n) =  m.sa;            // ...
        sm.block(a.n,a.n,b.n,b.n) = m.sb;         // ... 
        sm.block(0,a.n,a.n,b.n) = m.sab;          // ...
        sm.block(a.n,0,b.n,a.n) = m.sab.adjoint();

#ifdef DEBUG
	std::cout << "sm:" << std::endl;
        out_complex_matrix(sm,std::cout);
#endif
        
        ig = -sm;                                 // inverse-green function
        ig.block(0,0,a.n,a.n) += a.he - a.sl;     // ...
        ig.block(a.n,a.n,b.n,b.n) += b.he - b.sr; // ...
        
        dsm = sm;                                 // scatterer and lead self energy difference
        dsm.block(0,0,a.n,a.n) += -a.sr;          // ... 
        dsm.block(a.n,a.n,b.n,b.n) += -b.sl;      // ...

#ifdef DEBUG
	std::cout << "ig:" << std::endl;
        out_complex_matrix(ig,std::cout);
	std::cout << "dsm:" << std::endl;
        out_complex_matrix(dsm,std::cout);
#endif

        
        //g = ig.inverse();                       // calculating green function
        //k = g*dsm;                              // calculating lippman-swinger matrix
        k = ig.lu().solve(dsm);                    // more efficient..
        
#ifdef DEBUG
	std::cout << "k:" << std::endl;
        out_complex_matrix(k,std::cout);
#endif


        kaa = k.block(0,0,a.n,a.n);               // ... and its blocks
        kbb = k.block(a.n,a.n,b.n,b.n);           // ...
        kab = k.block(0,a.n,a.n,b.n);             // ...
        kba = k.block(a.n,0,b.n,a.n);             // ...

        raa = a.iwlo*kaa*a.wro;                   // calculating reflection and trnsmission matricies
        rbb = b.iwro*kbb*b.wlo;                   // ...
        tab = a.iwlo*kab*b.wlo;                   // ...
        tba = b.iwro*kba*a.wro;                   // ...

        ano = raa.rows();                         // dimensionality if relction matricies (open channels)
        bno = rbb.rows();                         // ...
        
        s.resize(ano+bno,ano+bno);                // constructiong scattering matrix
        s.block(0,0,ano,ano) = raa;               // ...
        s.block(ano,ano,bno,bno) = rbb;           // ...
        s.block(0,ano,ano,bno) = tab;             // ...
        s.block(ano,0,bno,ano) = tba;             // ...

        coefaa = raa.array().abs2().sum();        // reflection from lead-a to lead-a
        coefbb = rbb.array().abs2().sum();        // reflection from lead-b to lead-b
        coefab = tab.array().abs2().sum();        // transmission from lead-b to lead-a
        coefba = tba.array().abs2().sum();        // transmission from lead-a to lead-b

/*        
        ttab = tab*tab.adjoint();
        Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(ttab);  // solving eigenvalue problem
	ttabe = eigen_solver.eigenvalues();                                                                                    // eigenvalues
	ttabev = eigen_solver.eigenvectors();                                                                                    // eigenvalues
        
        xaar = a.wro + kaa*a.wro;                 // lead-a wf induced by right mov wf in lead-a
        xbar = kba*a.wro;                         // lead-b wf induced by right mov wf in lead-a

        xabl = kab*b.wlo;                         // lead-a wf induced by left mov wf in lead-b
        xbbl = b.wlo +kbb*b.wlo;                  // lead-b wf induced by left mov wf in lead-b
        
        xmr = - m.ihvae*xaar - m.ihvbe*xbar;      // scatterer wf induced by right mov wf in lead-a
        xml = - m.ihvae*xbbl - m.ihvbe*xabl;      // scatterer wf induced by left mov wf in lead-b
        
        yml = xml*ttabev;                         // scatterer wf for scatterer channels from lead-b to lead-a        
*/        
        maxcoefab = tab.array().abs2().maxCoeff();     // max transmission from lead-b to lead-a
        maxcoefba = tba.array().abs2().maxCoeff();     // max transmission from lead-a to lead-b

        
        id = s; id.setIdentity();                                    // calculate checknorm
        checknorm = ((s*s.adjoint()-id).array().abs()).maxCoeff();   // ... max abs (S*S'-I)

#ifdef DEBUG
        std::cout << "s*s':" << std::endl;
        out_ssd(e,s,std::cout);
#endif
        return true;
}


#ifndef PARALLEL
int main(int argc, char* argv[])
{
try
{
        lead a,b;
        scatterer m;
        
        double emax,emin,ef,tol;
        int n;
        bool sqrte,advanced;
        
        get_amb(argc,argv,a,m,b,emin,emax,ef,n,tol,sqrte,advanced);
        
        if (sqrte) { 
    	    emin = std::copysign(std::sqrt(std::abs(emin)),emin);
    	    emax = std::copysign(std::sqrt(std::abs(emax)),emax);
        }
	
	a.advanced = advanced;
	b.advanced = advanced;
	a.label = "a";
	b.label = "b";
        a.tol = tol;
        b.tol = tol;
        scattering s;
        double de = (emax-emin)/(n > 1 ? n : 0.5);
        
        
        out_coef(ef,s,std::cout,true);
        int try_catch_error_count = 0;
        for(double e = emin; e < emax+0.5*de; e+=de)
	{
            try
            {
    	    if (sqrte) { 
    		double ee = std::copysign(e*e,e);
    		s.transport(ee+ef,a,m,b); 
    		out_coef(ee,s,std::cout,false,advanced);
    	    }
    	    else {
    		s.transport(e+ef,a,m,b);
    		out_coef(e,s,std::cout,false,advanced);
    	    }
    	    
            
            } catch(const std::exception & e)
            {
                std::cerr << e.what() << std::endl;
                std::cerr << "Attempts to continue..." << std::endl;
                try_catch_error_count++;
                if(try_catch_error_count>10) throw "Too many errors!";
            }
        }
        
} catch(const char* message)
{
    std::cerr << message << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
catch(const std::exception & e)
{
    std::cerr << e.what() << std::endl;
    std::cerr << "(exit)" << std::endl;
    return EXIT_FAILURE;
}
    return EXIT_SUCCESS;
            

}
#else

#include "mpi.h"

void mpi_send_string(const std::string &str_msg, const int& to )
{
    const char* msg = str_msg.c_str();
    int len = str_msg.length();
    MPI::COMM_WORLD.Send( &len, 1,   MPI::INT,  to, ( int )1 );
    MPI::COMM_WORLD.Send(  msg, len, MPI::CHAR, to, ( int )1 );
}
            
void mpi_receive_string(std::string &str_msg, int& from, int buffer_size = 200 )
{
    int len;
    MPI::Status status;
    char* msg = new char[buffer_size];
    MPI::COMM_WORLD.Recv( &len, 1, MPI::INT,from,MPI::ANY_TAG,  status );
    MPI::COMM_WORLD.Recv(  msg, len, MPI::CHAR, from, MPI::ANY_TAG );
    msg[ len ] = '\0';
    str_msg = msg;
    delete [] msg;
}

int main(int argc, char* argv[])
{
try
{
        int mpi_id;
        int mpi_size;
        std::cerr << "MPI parallel execution will be performed." << std::endl;
        MPI::Init ( argc, argv );
        mpi_id = MPI::COMM_WORLD.Get_rank();
        mpi_size = MPI::COMM_WORLD.Get_size();
      
      
        lead a,b;
        scatterer m;
        
        double emax,emin,ef,tol;
        int n;
        bool advanced;
        
        get_amb(argc,argv,a,m,b,emin,emax,ef,n,tol,advanced);
        
        a.tol = tol;
        b.tol = tol;
        scattering s;
        double de = (emax-emin)/(n > 1 ? n : 0.5);

        int ne = n+1;


        if(ne < mpi_size) throw std::runtime_error("There are more nodes than energy points.");

        if(mpi_id==0)
        out_coef(ef,s,std::cout,true,advanced);
        for(int i=0; i < ne; i++)
        if(mpi_id == i%mpi_size)
        {
            double e = emin + double(i)*de;
            s.transport(e+ef,a,m,b);

            std::ostringstream os;
            out_coef(e,s,os,false,advanced);
            std::string line = os.str();
            
            if(mpi_id != 0)
                mpi_send_string(line,0);
            else
            {
                std::cout << line;
                int rest = mpi_size;
                if(n/mpi_size == i/mpi_size) rest = n%mpi_size+1;
                for(int id = 1; id < rest; id++)
                {
                    std::string received_line;
                    mpi_receive_string(received_line,id);
                    std::cout << received_line;
                }
            
            }
            
        }
         
        MPI::Finalize();
          
        
} catch(const char* message)
{
    std::cerr << message << std::endl;
    std::cerr << "(exit)" << std::endl;
    MPI::Finalize();
    return EXIT_FAILURE;
}
catch(const std::exception & e)
{
    std::cerr << e.what() << std::endl;
    std::cerr << "(exit)" << std::endl;
    MPI::Finalize();
    return EXIT_FAILURE;
}
    return EXIT_SUCCESS;
}
#endif

